#include <R.h>
#include "Rmath.h"

extern "C" {

// Algorithm AS 204 Appl. Statist. (1984) Vol. 33, No.3
// ruben evaluates the probability that a positive definite quadratic form in Normal variates is less than a given value
  
  void ruben(double *lambda, int *mult, double *delta, int *n, double *c, double *mode, int *maxit, double *eps, double *dnsty, int *ifault, double *res) {
    
    int i,k,m,j;
    double pnorm(double q, double mean, double sd, int lower_tail, int log_p);
    double ao, aoinv, z, bbeta, eps2, hold, hold2, sum, sum1, dans, lans, pans, prbty, tol;
    double *gamma, *theta, *a, *b;
    gamma = new double[n[0]];
    theta = new double[n[0]];
    a = new double[maxit[0]];
    b = new double[maxit[0]];
    
    
    
    if ((n[0]<1) || (c[0]<=0) || (maxit[0] <1) || (eps[0]<=0.0)) {
      res[0] = -2.0;
      ifault[0] = 2;
      delete[] gamma;
      delete[] theta;
      delete[] a;
      delete[] b;
      return;
    } else {
      tol = -200.0;
    
      // Preliminaries
      sum = lambda[0];
      bbeta = sum;
      
      for (i=1;i<=n[0];i++) {
	hold = lambda[i-1];
	if ((hold<=0.0) || (mult[i-1]<1) || (delta[i-1]<0.0)) {
	  res[0] = -7.0;
	  ifault[0] = -i;
	  delete[] gamma;
	  delete[] theta;
	  delete[] a;
	  delete[] b;
	  return;
	}	
	if (bbeta > hold) bbeta = hold; // calcul du max des lambdas
	if (sum < hold) sum = hold;    // calcul du min des lambdas
      }
      
  
    if (mode[0] > 0.0) {
      // if ((2.0/(1.0/bbeta+1.0/sum))>1.8*sum) bbeta = sum; // comme dans NAG : methode avec betaA
            bbeta = mode[0]*bbeta;
    } else {
      bbeta = 2.0/(1.0/bbeta+1.0/sum);  // methode avec betaB
    }

    k = 0;
    sum = 1.0;
    sum1 = 0.0;
    for (i=1;i<=n[0];i++) {
      hold = bbeta/lambda[i-1];
      gamma[i-1] = 1.0 - hold;
      sum = sum*R_pow(hold,mult[i-1]); //???? pas sur ..
      sum1 = sum1 + delta[i-1];
      k = k + mult[i-1];
      theta[i-1] = 1.0;
    }
    
    ao = exp(0.5*(log(sum)-sum1));
    if (ao <= 0.0) {
      res[0] = 0.0;
      dnsty[0] = 0.0;
      ifault[0] = 1;
    } else { // evaluate probability and density of chi-squared on k degrees of freedom. The constant 0.22579135264473 is ln(sqrt(pi/2))
      z = c[0]/bbeta;
      
      if ((k%2) == 0) { // k est un entier donc on regarde si k est divisible par 2: k == (k/2)*k 
	i = 2;
	lans = -0.5*z;
	dans = exp(lans);
	pans = 1.0 - dans;
      } else {
	i = 1;
	lans = -0.5*(z+log(z)) - 0.22579135264473;
	dans = exp(lans);
	pans = pnorm(sqrt(z),0.0,1.0,1,0) - pnorm(-sqrt(z),0.0,1.0,1,0); 
      }
      
      k = k-2;
      for (j=i;j<=k;j=j+2) {
	if (lans < tol) {
	  lans = lans + log(z/(double)j);
	  dans = exp(lans);
	} else {
	  dans = dans*z/(double)j;
	}
	pans = pans -dans;
      }
      
      // evaluate successive terms of expansion
      
      prbty = pans;
      dnsty[0] = dans;
      eps2 = eps[0]/ao;
      aoinv = 1.0/ao;
      sum = aoinv - 1.0;
    

    for (m=1;m<=maxit[0];m++) {
      sum1 = 0.0;
      for (i=1;i<=n[0];i++) {
	hold = theta[i-1];
	hold2 = hold*gamma[i-1];
	theta[i-1] = hold2;
	sum1 = sum1 + hold2*mult[i-1]+m*delta[i-1]*(hold-hold2);
      }
      sum1 = 0.5*sum1;
      b[m-1] = sum1;
      for (i=m-1;i>=1;i--) {
	sum1 = sum1 + b[i-1]*a[m-i-1]; 
      }
      sum1 = sum1/(double)m;
      a[m-1] = sum1;
      k = k + 2;
      if (lans < tol) {
	lans = lans + log(z/(double)k);
	dans = exp(lans);
      } else {
	dans = dans*z/(double)k;
      }
      pans = pans - dans;
      sum = sum - sum1;
      dnsty[0] = dnsty[0] + dans*sum1;
      sum1 = pans*sum1;
      prbty = prbty + sum1;
      if (prbty<(-aoinv)) {
	res[0] = -3.0;
	ifault[0] = 3;
	return;
      }
      if (fabs(pans*sum) < eps2) {
	if (fabs(sum1) < eps2) {
	  ifault[0] = 0;
	  
	  m = maxit[0]+1;
	  break;

	}
      }
    }

    ifault[0] = 4;
    dnsty[0] = ao*dnsty[0]/(bbeta+bbeta);
    prbty = ao*prbty;
    if (prbty<0.0 || prbty>1.0) {ifault[0] = ifault[0] + 5;
    } else {
      if (dnsty[0]<0.0) ifault[0] = ifault[0] + 6;
    }
    res[0] = prbty;
    }

    delete[] gamma;
    delete[] theta;
    delete[] a;
    delete[] b;
    return;
    }
    
  }
}

